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RANDOM SAMPLING FOR MULTIVARIATE BERNOULLI VARL\BLES 



BACKGROUND 

1 . Field of the Invention 

[0001] This invention relates to networks. In particular, the invention relates to network 
reliability. 

2. Descriptionof Related Art 

[0002] In certain studies of reliability of computer networks, it is useful to be able to 
generate samples of a Bernoulli random vector whose components model the current state 
of some network element. For example, the individual Bernoulli variables may correspond 
to the links in the network, with a value of 1 indicating the link is functional or "up", and a 
value of 0 indicating it is non-functional or "down". 

[0003] Existing work on network reliability usually assumes that the Bernoulli variables 
corresponding to the links are independent. This assumption means that the random 
samples can be obtained independently for the individual components. Under this 
assumption, if the probability of each link being up or down is specified, then the sampling 
problem is trivial. However, the independence assumption is not realistic for certain types 
of networks, such as those using free-space optical (FSO) links. In these networks, there 
are common factors that affect the states of the links. For example, weather conditions 
may affect links within a locality. If a particular FSO link is down due to heavy fog, then 
it is likely that nearby links are down as well. The states of the links in these types of 
networks are, therefore, correlated and/or dependent on each other. To characterize this 
correlation, a covariance matrix is used in addition to the means of the link variables. The 
necessity of providing a covariance matrix for a Bernoulli variables makes it difficult to 
generate random samples of the network state. 

[0004] In addition, the computational requirements for the joint probability density of such 
random vector are enormous. This is because such a joint density is a discrete density 
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having 2" components where n is the length of the random vector. For a vector of even 
moderate length, the number of components making up the discrete density is extremely 
large, over 1 billion for n = 30. 

[0005] Therefore, there is a need to have an efficient technique to generate a random 
vector for multivariate Bernoulli variables. 
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BRIEF DESCRIFnON OF THE DRAWINGS 

[0006] The features and advantages of the present invention will become apparent from the 
following detailed description of the present invention in which: 

[0007] Figure 1 is a diagram illustrating a system in which one embodiment of the 

invention can be practiced. 

[0008] Figure 2 is a diagram illustrating a simulator shown in Figure 1 according to one 
embodiment of the invention. 

[0009] Figure 3 is a diagram illustrating a random sampler shown in Figure 2 according to 

one embodiment of the invention. 

[0010] Figure 4 is a diagram illustrating a co variance generator shown in Figure 3 
according to one embodiment of the invention. 

[0011] Figure 5 is a flowchart illustrating a process to generate samples for a Bernoulli 
distribution according to one embodiment of the present invention. 
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DESCRIPTION 

[0012] One embodiment of the present invention generates random samples for a Bernoulli 
distribution. A first covariance matrix is generated using a desired mean vector and a 

desired covariance matrix of the Bernoulli distribution. A normal vector is constructed 
using the desired mean vector and the first covariance matrix. A sampling v^tor is 
generated using the normal vector and a threshold vector. The sampling vector has the 
desired mean vector and the desired covariance matrix. 

[0013] In the following description, for purposes of explanation, numerous details are set 
forth in order to provide a thorough understanding of the present invention. However, it 
will be apparent to one skilled in the art that these specific details are not required in order 
to practice the present invention. In other instances, well-known electrical structures and 
circuits are shown in block diagram form in order not to obscure the present invention. 

[0014] Figure 1 is a diagram illustrating a computer system 100 in which one embodiment 
of the invention can be practiced. The computer system 100 includes a processor 110, a 
host bus 120, a memory control hub (MCH) 130, a system memory 140, an input/output 
control hub (ICH) 150, a mass storage device 170, and input/output devices 180i to 180k. 

[0015] The processor 110 represents a central processing unit of any type of architecture, 
such as embedded processors, micro-controllers, digital signal processors, superscalar 
computers, vector processors, single instruction multiple data (SIMD) computers, complex 
instruction set computers (CISC), reduced instruction set computers (RISC), very long 
instruction word (VLIW), or hybrid architecture. The host bus 120 provides interface 
signals to allow the processor 1 10 to communicate with other processors or devices, e.g., 
the MCH 130. The host bus 120 may support a imi-processor or multiprocessor 
configuration. The host bus 120 may be parallel, sequential, pipelined, asynchronous, 
synchronous, or any combination thereof. 

[0016] The MCH 130 provides control and configuration of memory and input/output 
devices such as the system memory 140 and the ICH 150. The MCH 130 may be 
integrated into a chipset tiiat integrates multiple fimctionalities such as the isolated 
execution mode, host-to-peripheral bus interface, memory control. For clarity, not all the 
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peripheral buses are shown. It is contemplated that the system 100 may also include 
peripheral buses such as Peripheral Component Interconnect (PCI), accelerated graphics 
port (AGP), Industry Standard Architecture (ISA) bus, and Universal Serial Bus (USB), 
etc. 

[0017] The system memory 140 stores system code and data. The system memory 140 is 
typically implemented with dynamic random access memory (DRAM) or static random 
access memory (SRAM). The system memory may include program code or code 
segments implementing one embodiment of the invention. The system memory includes a 
network simulator 145. The network simulator 145 may also be implemented by 
hardware, software, firmware, microcode, or any combination thereof. The system 
memory 140 may also include other programs or data which are not shown, such as an 
operating system. 

[0018] The ICH 150 has a number of fiuictionalities that are designed to support I/O 
functions. The ICH 150 may also be integrated into a chipset together or separate from the 
MCH 130 to perform I/O functions. The ICH 150 may include a number of interface and 
I/O functions such as PCI bus interface, processor interface, interrupt controller, direct 
memory access (DMA) controller, power management logic, timer, universal serial bus 
(USB) interface, mass storage interface, low pin count (LPC) interface, etc. 

[0019] The mass storage device 170 stores archive information such as code, programs, 
files, data, applications, and operating systems. The mass storage device 170 may include 
compact disk (CD) ROM 172, floppy diskettes 174, and hard drive 176, and any other 
magnetic or optic storage devices. The mass storage device 170 provides a mechanism to 
read machine-readable media. 

[0020] The I/O devices 180i to 1 80k may include any FO devices to perform VO 
functions. Examples of I/O devices ISOi to 180k include controller for input devices (e.g., 
keyboard, mouse, trackball, pointing device), media card (e.g., audio, video, graphics), 
network card, and any other peripheral controllers. 

[0021] The present invention may be implemented by hardware, software, firmware, 
microcode, or any combination thereof. When implemented in software, firmware, or 
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microcode, the elements of the present invention are the program code or code segments to 
perform the necessary tasks. A code segment may represent a procedure, a function, a 
subprogram, a program, a routine, a subroutine, a module, a software package, a class, or 
any combination of instructions, data structures, or program statements. A code segment 
may be coupled to another code segment or a hardware circuit by passing and/ or receiving 
information, data, arguments, parameters, or memory contents. Information, argxmients, 
parameters, data, etc. may be passed, forwarded, or transmitted via any suitable means 
including memory sharing, message passing, token passing, network transmissionj etc. The 
program or code segments may be stored in a processor readable medium or transmitted by 
a computer data signal embodied in a carrier wave, or a signal modulated by a carrier, over 
a transmission medium. The "processor readable medium" may include any medium that 
can store or transfer information. Examples of the processor readable medium include an 
electronic circuit, a semiconductor memory device, a ROM, a flash memory, an erasable 
ROM (EROM), a floppy diskette, a compact disk CD-ROM, an optical disk, a hard disk, a 
fiber optic medium, a radio frequency (RF) link, etc. The computer data signal may 
include any signal that can propagate over a transmission medium such as electronic 
network channels, optical fibers, air, electromagnetic, RF links, etc. The code segments 
may be downloaded via computer networks such as the Internet, Intranet, etc. 

[0022] It is noted that the invention may be described as a process which is usually 
depicted as a flowchart, a flow diagram, a structure diagram, or a block diagram. Although 
a flowchart may describe the operations as a sequential process, many of the operations 
can be performed in parallel or concurrently. In addition, the order of the operations may 
be re-arranged. A process is terminated when its operations are completed. A process may 
correspond to a method, a function, a procedure, a subroutine, a subprogram, etc. When a 
process corresponds to a function, its termination corresponds to a return of the function to 
the calling function or the main function. 

[0023] Figure 2 is a diagram illustrating the simulator 145 shown in Figure 1 according to 
one embodiment of the invention. The simulator 145 includes a network modeler 210, a 
random sampler 220, a reliability modeler 230, and a network performance evaluator 240. 
The simulator 145 may be implemented by software, hardware, or any combination 
thereof. Each of the network modeler 210, the random sampler 220, the reliability modeler 
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230, and the network performance evaluator 240 may be a software module, a function, a 
subprogram, a method, or a hardware element, or any combination of hardware and 
software. 

[0024] The network modeler 210 creates a model for the network whose rehability is to be 
studied. The network modeler 210 receives the network parameters such as number of 
hubs, number of links, weather conditions, etc. The network model generated by the 
network modeler 210 may be any type of network and network configuration. In one 
embodiment, the network is a FSO network. An example of this network is a FSO 
network 250. The network 250 may be located in a metropolitan area. In the FSO network 
250, hubs are connected via fiberless optical links. At a transmitter at a hub, electrical 
signal representing the data or information is converted into a beam of light by a laser or 
optical transmitter. The light is then sent through a lens that increases its divergence. This 
wide beam of Ught is then sent through the air to the receiving unit at the other hub. As the 
Ught enters the receiver, it goes through several lenses and is focused on a small detector. 
The detector then converts the Ught back into electrical signals for the network equipment 
to use. 

[0025] hi the network 250, there are N hubs 260i to 260n. A hub is connected to another 
hub through an optical link. A hub may have one or more transmitters and/or receivers so 
that it can transmit to or receive from multiple hubs. A hub 260k is connected to a hub 260j 
via a link 270kj. For example, hub 260i and 26O2 is connected by optical link 270i2. There 
are M optical links 270i to 270m. Since these links exist within a locality of a metropolitan 
area through the air medium, their operational states are highly correlated and dependent. 
For example, in a sunny day, the optical links are highly operational and their states may 
be characterized as operational or "up". On the other hand, in a foggy day, these links may 
become non-operational and then- states may be characterized as non-operational or 
"down". 

[0026] The random sampler 220 receives the model parameters from the network modeler 
210 to generate random samples to the reliability modeler 230. The model parameters may 
include statistical parameters of the states of the optical links m the network 250. The 
statistical parameters may include the distribution type, the desired mean, and the desired 
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covaiiance matrix of the distribution. Since the state of an optical Unk can assume one of 
two values: up or down, the distribution may be of a Bernoulli distribution. 

[0027] The reliability modeler 230 creates a reliability model to study the reliability of the 
network 250 as modeled by the network modeler 210. The reliability model may be any 
convenient model including combinatorial models, Markov models, or any other 
probabiUstic models. The reliability model 230 receives the random samples from the 
random sampler 220 and the network information from the network modeler 210 to apply 
to the reliability model. 

[0028] The network performance evaluator 240 evaluates the network performance based 
on the reliability as provided by the reliability modeler. The network performance 
evaluator 240 may provide useful information regarding network behavior so that 
appropriate network configuration may be employed for a wide range of network 
conditions. 

10029] The random sampler 220, therefore, is an important component of the simulator 
145 to provide realistic information regarding the conditions of the network 250 so that its 
reUability may be accurately predicted or evaluated. Since the states of the optical links 
are correlated and/or dependent, the construction of random samples based on a Bernoulli 
distribution cannot be performed on an element by element basis. 

[0030] The random sampler 220 generates the samples to have the desired mean vector \x 
and the desired covariance matrix S. The procedure to do this is explained in the 
following. First, the scalar case is considered. The scalar case provides a formula to 
compute tiie threshold elements in a threshold vector to be used in the generation of the 
random samples. Then, the multivariate case is considered. 

The scalar case: 

[0031] Let X be a Bernoulli random variable witii mean n and variance cr^. Then P(X=1) 
^ [I, and = |J.(l-li). Given such an X, let Y be a Gaussian random variable with the same 
mean and variance. Thus, Y ~ N(n,a^) where N denote the normal distribution. If a 
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threshold x is chosen, a new Bernoulli random variable can be defined by the following 
rule: 

Z = liff (if and only if) Y > x, Z = 0 otherwise ( 1 ) 

[0032] It is desired to choose x such that Z has a mean jj, and variance a^. Since cp^ = 

\x), it is enough to choose x to get the desired mean. Note that Z has mean |j, iff P(Z=1) = p, 

which, in turn, is equivalent to P(Y> x) = ]x. But: 



P(Y>x) = P((Y-|Li)/a) > (x-^i)/a) =. -L I e 2 dt (2) 



[0033] Therefore, it is necessary to find x so that the integral on the right half side (RHS) is 
\i. To find the solution, recall that the error function is defined as: 

2 z 2 

erf(z)= -^f e-^ dt (3) 

[0034] From this, the expression in (2) may be written as: 

_l^]e-d. =io-erf(i)) (4) 

[0035] Using equation (4) and equate (2) to jti, it follows: 

i(,_erf(ll|)) = , (5) 

[0036] Therefore, the solution is: 

X = )a + a V2 inverf( 1 -2\x) (6) 
The multivariate case: 
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[00371 Consider a random vector of Bernoulli variables X = (Xi, . . Xn) having a mean 
vector = (ni, . . M and covariance matrix E = (Eij), where N is a positive integer 
denoting the number of links in the network. Note that = \Xi(l-\ii). A multivariate 
Gaussian distribution is constructed from the mean vector (a and the covariance matrix S. 
Let Y = (Yi, . . Yn) be a random vector drawn from this Gaussian or normal distribution. 

Using equation (6), let t = (xi Xn) be the threshold vector. A new vector of Bernoulli 

variables, Z = (Zi, . . Zn) is defined by the following rule: 

Zi = 1 iff Yi> Xi, Zi= 0 otherwise (7) 

[00381 Each Bernoulli variable Zi has the same mean and variance as the corresponding X 
and Yi. However, generally the covariance matrix of Z is not E. The diagonal elements 
will be those of S because Xj, Yi, and Zi have the same variance, but the non-diagonal 
elements will not be the same. To construct Z to have the desired covariance matrix, 
another covariance matrix will be generated. This covariance matrix will be used as the 
basis for the construction of a second normal distribution from which the thresholding rule 
of (7) can be applied. The procedure to generate this covariance matrix will be explained 
later. 

[00391 For now, assuming that this can be done. Observe that the random samples of Z 
can be obtained by first obtaining samples of Y and then computing the correspondmg 
values of Z. Originally, only the mean vector and the covariance matrix for Z are known, 
and Z cannot be sampled because the joint densities are not known. By relating to the 
Gaussian vector Y, the joint density for Z can be known, namely: 

fv..,ZN(l, 0,...,1) =P(Zi-l, Z2=0,...,Zn=1) 

= P(Yi>xi , Y2 < X2, . . . , Yn > Xn), etc. (8) 

[00401 It is well known that the usual multivariate normal distribution corresponding to \i 
and 2 is the maximum entropy distribution having the given mean and the given 
covariance matrix. 

Constructing another covariance matrix: 
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[0041] Instead of using the desired covariance matrix S, the procedure starts with 
construction of a matrix S, chosen in such a way that the covariance matrix of the random 
vector Z is the desired covariance matrix S. In essence, the procedure is as follows: 

[0042] Step 1: Write down Ihe expression for the covariance of Zi and Zj and set it equal 
to Sij. The expression for the covariance of Zi and Zj depends on Hi, ^ij, and on the non- 
diagonal elements Sij of the matrix S. Note that the diagonal elements sa = ^(1-^0- 

[0043] Step 2: Solve the equation numerically for sy. 

10044] Step 3: Repeat steps 1 and 2 for all all elements of the covariance matrix S 
can then be obtained. 

[0045] To write the equation for sy, note that: 
cov(Zi,Zj) = E(ZiZj)-jiiHj 

= P(Zi=landZj=l)-KiiM.j 

= P(Yi>Xi and Yj>Tj) - mm (9) 
[0046] It can be written: 

P(Yi>TiandYj>Xj) = 

o^J: 2(1 -p2) "^i^j 

[0047] where a = l^ra^aj ^l-p^ 

p=Sij/cTiaj (lib) 
[0048] Introducing changes of variables: 

u = (x - ^ii)/ai and V = (y - Hj)/cj (12) 
[0049] Then the RHS of equation (10) becomes: 



11 



^ f J exp( ^— ^ - 2 ffuv + )rfMifv 



(13) 



10050] Let c > 0 be defined by 

c2 = 2(l-p^) 
[00511 and define: 

Si = (Ti - )Ji)/ai for all i = 1, N 
[0052] Then introduce new coordinates p and q by: 

p = (x+y)/c V2 and q = (y-x)/c V2 

[0053] and finally define variables P and Q by: 

' r- andQ= ^ ^' 

cV2 cV2 

[0054] The integral expression in equation (10) becomes: 



/l _ 2 « P+Q-P 
^ p+Q+P 



(14) 



(15) 



(16) 



(17) 



(18) 



[0055] The inner integral can be expressed in terms of the error function. The integral 
expression finally becomes: 



F(p)=iiP fe-P (l-P)(erf(Vr?p(Q + P + P))-erf(^^(Q-p + P)))dp 

2V7C ^ 

(19) 



[0056] The equation to be solved therefore is: 
F(p) = 2ij + HiHj 



(20) 
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[00571 The eqviation (20) can be solved numerically, such as a modified secant method. 
From values of p, it is straightforward to obtain s = p <t, (Tj as shown in equation (lib). 

[00581 Figure 3 is a diagram illustrating the random sampler 220 shown in Figure 2 
according to one embodiment of the invention. The random sampler 220 includes a 
covariance generator 310, a normal vector generator 320, and a thresholder 330. 

[0059] The covariance generator 310 receives a desired mean vector \x = (m, . . ., Hn) and a 
desired covariance matrix E = (Sy) where i = 1, N andj = 1, . . N. N is a positive 
integer and indicates the number of optical links in the network. The desired mean vector 
|4, and the desired covariance matrix S are the mean vector and covariance matrix of a 
multivariate Bernoulli distribution. Typically, the desired mean vector n and the desired 
covariance matrix Z are specified as part of the simulation to study the reliability of the 
network. The covariance generator 310 generates a covariance matrix S and a threshold 
vector T = (ti, . . ., Xn). 

[00601 The normal vector generator 320 generates a normal vector Y = (Yi, . . Yn) based 
on the desired mean vector// and the covariance matrix S. The normal vector Y has a 
multivariate Gaussian distribution having the mean vector equal to the desired mean vector 
and the covariance matrix equal to S. 

[00611 The thresholder 330 applies the threshold vector t to generate the desired samples Z 
based on the normal vector Y according to the threshold rule shown in equation (7). The 
desired samples z are samples of a Bernoulli distribution having a mean// and a covariance 

matrix E. 

[0062] Figure 4 is a diagram illustrating the covariance generator 310 shown in Figure 3 
according to one embodiment of the invention. The covariance generator 310 includes a 
normal vector generator 410, a threshold vector calculator 420, an integral expression 
generator 430, a non-diagonal element generator 440, and a diagonal element generator 
450. 
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[0063] The normal vector generator 410 generates a Gaussian distribution from the desired 
mean vector |j, and the desired covariance matrix 2. The normal vector generator 410 has 
N normal elements. Each of the normal elements at a vector index k having a mean Hk and 
a variance . The mean is equal to a desired mean element of the desired mean vector 
(i at the vector index k. 

[0064] The threshold vector calculator 420 calculates the threshold elements of the 
threshold vector x = (xi, . . xn) as follows: 

Ti = Hi + Oi >/2 inverf(l-2|Xi) (21) 

[0065] where inverf is an inverse error function as shown in equation (6). 

[0066] The integral expression generator 430 generates an integral expression F(p) using 
the threshold vector t. The integral expression generator 430 forms the variables p, c, P, 
andQ as given in equations (1 lb), (14), and (17). 

[0067] Note that P is the lower integral limit of the integral expression F(p). The upper 
integral limit is infinity. Then, the integral expression generator 430 generates the integral 

expression F(p) as shown in equation (19). 

[0068] The non-diagonal element generator 440 obtains the first non-diagonal element Sij 
using the integral expression F(p), the means m and \Xj, and a desired non-diagonal element 
Zij of the desired covariance matrix. The non-diagonal element generator 440 includes a 
RHS generator 442 and an equation solver 444. 

[00iS9] The RHS generator 442 generates the RHS of the equation for the integral 
expression provided by the integral expression generator 430. The RHS generator 442 
receives the desired mean vector \i and the desired covariance matrix S and calculates the 
RHS quantity gij as: 
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[00701 The equation solver 444 equates the integral expression F(p) with the RHS quantity 
gij to form the equation: 

F(p) = gij (23) 

[0071] In this equation, for given values of Xi, xj, ^i, |J.j, Oi, and aj, the quantity p can be 
determined. Once p is determined, the non-diagonal element Sij can be derived directly. 
The integral equation (23) can be solved numerically by a number of techniques. In one 
embodiment, the integral equation (23) is solved by using a modified secant method as 
implemented in Mathematica's FindRoot fimction, 

100721 The diagonal element generator 450 computes the diagonal elements Sjj of the 
covariance matrix S by equating sjj to the desired diagonal element 2jj = /Jj (1 - //^ ) of the 
desired covariance matrix S. The diagonal elements and the non-diagonal elements define 
the covariance matrix to be provided to the normal vector generator 320 as shown in 
Figure 3. 

[00731 Figure 5 is a flowchart illustrating a process 500 to generate samples for a 
Bernoulli distribution according to one embodiment of the present invention. 

[0074] Upon START, the process 500 initializes the row and column indices (ij) for the 
covariance matrix (Block 520). Next, the process 500 generates the threshold vector x = 
(xi, . . ., Xn) using equation (21) (Block 530). Next, the process 500 generates the integral 
expression F(p) using the equation (19) (Block 540). Then, the process 500 generates the 
RHS quantity gy using the equation (22) (Block 545). Next, the process 500 forms the 
integral equation and solve the integral equation for p and obtains the elements Sjj for the 
first covariance matrix (Block 550). In addition, the process 500 computes the diagonal 
elements sh of the first covariance matrix. 

[00751 Then, the process 500 determines if all elements at all row and column indices have 
been processed (Block 560). If not, the process 500 updates the row and/or column indices 
(Block 565) and goes to Block 540. Otherwise, the process 500 generates a normal vector 
Y = (Yi, . . Yn) fi"om the desired mean vector and the first covariance matrix S (Block 
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570). Next, the process 500 calculates the sample Z by thresholding using the threshold 
vector T (Block 580) and is then terminated. 

[0076] While this invention has been described with reference to illustrative embodiments, 
this description is not intended to be construed in a limiting sense. Various modifications 
of the illustrative embodiments, as well as other embodiments of the invention, which are 
apparent to persons skilled in the art to which the invention pertains are deemed to lie 
within the spirit and scope of the invention. 
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